Load the necessary libraries
library(car) # for regression diagnostics
library(broom) # for tidy output
library(broom.mixed) # for tidy output
library(ggfortify) # for model diagnostics
library(sjPlot) # for outputs
library(knitr) # for kable
library(effects) # for partial effects plots
library(ggeffects) # for effects plots in ggplotjk
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(MuMIn) # for AICc
library(tidyverse) # for data wrangling
library(DHARMa) # for assessing dispersion etc
library(glmmTMB) # for glmmTMB
library(performance) # for diagnostic plots
library(see) # for diagnostic plots
library(patchwork) # for multiple plots
fanworms
In an attempt to understand the effects on marine animals of short-term exposure to toxic substances, such as might occur following a spill, or a major increase in storm water flows, a it was decided to examine the toxicant in question, Copper, as part of a field experiment in Hong Kong. The experiment consisted of small sources of Cu (small, hemispherical plaster blocks, impregnated with copper), which released the metal into sea water over 4 or 5 days. The organism whose response to Cu was being measured was a small, polychaete worm, Hydroides, that attaches to hard surfaces in the sea, and is one of the first species to colonise any surface that is submerged. The biological questions focused on whether the timing of exposure to Cu affects the overall abundance of these worms. The time period of interest was the first or second week after a surface being available.
The experimental setup consisted of sheets of black perspex (settlement plates), which provided good surfaces for these worms. Each plate had a plaster block bolted to its centre, and the dissolving block would create a gradient of [Cu] across the plate. Over the two weeks of the experiment, a given plate would have plain plaster blocks (Control) or a block containing copper in the first week, followed by a plain block, or a plain block in the first week, followed by a dose of copper in the second week. After two weeks in the water, plates were removed and counted back in the laboratory. Without a clear idea of how sensitive these worms are to copper, an effect of the treatments might show up as an overall difference in the density of worms across a plate, or it could show up as a gradient in abundance across the plate, with a different gradient in different treatments. Therefore, on each plate, the density of worms (#/cm2) was recorded at each of four distances from the centre of the plate.
Format of copper.csv data file
| COPPER | PLATE | DIST | WORMS | AREA | COUNT |
|---|---|---|---|---|---|
| .. | .. | .. | .. | .. | .. |
| COPPER | Categorical listing of the copper treatment (control = no copper applied, week 2 = copper treatment applied in second week and week 1= copper treatment applied in first week) applied to whole plates. Factor A (between plot factor). |
| PLATE | Substrate provided for polychaete worm colonisation on which copper treatment applied. These are the plots (Factor B). Numbers in this column represent numerical labels given to each plate. |
| DIST | Categorical listing for the four concentric distances from the center of the plate (source of copper treatment) with 1 being the closest and 4 the furthest. Factor C (within plot factor) |
| WORMS | Density (#/cm2) of worms measured. Response variable. |
copper <- read_csv("../public/data/copper.csv", trim_ws = TRUE)
glimpse(copper)
## Rows: 60
## Columns: 4
## $ COPPER <chr> "control", "control", "control", "control", "control", "control…
## $ PLATE <dbl> 200, 200, 200, 200, 39, 39, 39, 39, 34, 34, 34, 34, 36, 36, 36,…
## $ DIST <dbl> 4, 3, 2, 1, 4, 3, 2, 1, 4, 3, 2, 1, 4, 3, 2, 1, 4, 3, 2, 1, 4, …
## $ WORMS <dbl> 11.50, 13.00, 13.50, 12.00, 17.75, 13.75, 12.75, 12.50, 11.50, …
Let start by declaring the categorical variables and random effect as factors.
copper <- copper %>%
mutate(
COPPER = factor(COPPER),
PLATE = factor(PLATE),
DIST = factor(DIST)
)
In the current example, the response is worm density. Density (count per area) is a particularly hard variable to model against because it has both count and measurement properties. In some ways, it is easier to model the count against a Poisson (or Negative Binomial) distribution and then use an offset for Area so as to effectively model density. Unfortunately, we only have the densities (not the original counts and areas).
We therefore have six options (if density does not follow a Gaussian distribution - which it is unlikely to do):
Model formula: \[ y_i \sim{} \mathcal{Gamma}(\lambda_i, \phi)\\ ln(\lambda_i) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of copper, distance and their interaction on the number of number of worms. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual plates. \(\phi\) represents the shape parameter of the Gamma distribution.
Let start by exploring the patterns in worm density against both distance and the copper treatment.
ggplot(copper, aes(y = WORMS, x = DIST, fill = COPPER)) +
geom_boxplot()
Conclusions:
ggplot(copper, aes(y = WORMS, x = DIST, fill = COPPER)) +
geom_boxplot() +
scale_y_continuous(trans = scales::pseudo_log_trans())
Conclusions:
It might also be useful to see the raw data
g1 <- ggplot(copper, aes(y = WORMS, x = DIST, colour = COPPER)) +
geom_point(position = position_jitter(height = 0, width = 0.1)) +
scale_y_continuous()
g2 <- ggplot(copper, aes(y = WORMS, x = DIST, colour = COPPER)) +
geom_point(position = position_jitter(height = 0, width = 0.1)) +
scale_y_continuous(trans = scales::pseudo_log_trans())
g3 <- ggplot(copper, aes(y = WORMS, x = DIST, colour = COPPER)) +
geom_point(position = position_jitter(height = 0, width = 0.1)) +
scale_y_sqrt()
g1 + g2 + g3
Conclusions:
We can also explore the distance effects for each PLATE
ggplot(copper, aes(y = WORMS, x = as.numeric(PLATE), color = DIST)) +
geom_line() +
facet_wrap(~COPPER)
Conclusions:
Lets make the minimum half the next smallest value
copper %>%
filter(WORMS > 0) %>%
summarize(min = min(WORMS) / 2)
copper <- copper %>% mutate(WORMSp = ifelse(WORMS == 0, 0.125, WORMS))
Exploratory data analysis suggested that there might be an issue with homogeneity of variance above what might be expected to be addressed by a log-link. One way of dealing with this is to relax the homogeneity of variance assumption and instead indicate that variance (dispersion) should be calculated separately per group. This is supported in glmmTMB, however there has been an issue with calculating confidence intervals when this is the case. To address this (as of Sept 2020), the glmmTMB development team suggest installing glmmTMB from GitHub to get the latest fixes.
remotes::install_github("glmmTMB/glmmTMB/glmmTMB")
Whilst it might have been a bit difficult to determine whether or not normality was satisfied, it was more clear that variances were not equal between the treatments. glmmTMB allows us to fit models in which the variance is modelled against predictors.
copper.glmmTMB1a <- glmmTMB(WORMS ~ COPPER * DIST + (1 | PLATE),
dispformula = ~ COPPER * DIST,
data = copper,
family = gaussian(),
REML = TRUE,
control = glmmTMBControl(
optimizer = optim,
optArgs = list(method = "BFGS")
)
)
copper.glmmTMB1b <- glmmTMB(WORMS ~ COPPER * DIST + (DIST | PLATE),
dispformula = ~ COPPER * DIST,
data = copper,
family = gaussian(),
REML = TRUE,
control = glmmTMBControl(
optimizer = optim,
optArgs = list(method = "BFGS")
)
)
AICc(copper.glmmTMB1a, copper.glmmTMB1b)
## copper.glmmTMB1 <- update(copper.glmmTMB1, REML=TRUE)
Conclusions:
copper.glmmTMB2a <- glmmTMB(sqrt(WORMS) ~ COPPER * DIST + (1 | PLATE),
dispformula = ~ COPPER * DIST,
data = copper,
family = gaussian(),
REML = TRUE
)
copper.glmmTMB2b <- glmmTMB(sqrt(WORMS) ~ COPPER * DIST + (DIST | PLATE),
dispformula = ~ COPPER * DIST,
data = copper,
family = gaussian(),
REML = TRUE,
control = glmmTMBControl(
optimizer = optim,
optArgs = list(method = "BFGS")
)
)
AICc(copper.glmmTMB2a, copper.glmmTMB2b)
Conclusions:
copper.glmmTMB3a <- glmmTMB(WORMSp ~ COPPER * DIST + (1 | PLATE),
dispformula = ~ COPPER * DIST,
data = copper,
family = gaussian(link = "log"),
REML = TRUE
)
copper.glmmTMB3b <- glmmTMB(WORMSp ~ COPPER * DIST + (DIST | PLATE),
dispformula = ~ COPPER * DIST,
data = copper,
family = gaussian(link = "log"),
REML = TRUE,
control = glmmTMBControl(
optimizer = optim,
optArgs = list(method = "BFGS")
)
)
AICc(copper.glmmTMB3a, copper.glmmTMB3b)
Conclusions:
copper.glmmTMB4a <- glmmTMB(WORMSp ~ COPPER * DIST + (1 | PLATE),
dispformula = ~ COPPER * DIST,
data = copper,
family = Gamma(link = "log"),
REML = TRUE
)
copper.glmmTMB4b <- glmmTMB(WORMSp ~ COPPER * DIST + (DIST | PLATE),
dispformula = ~ COPPER * DIST,
data = copper,
family = Gamma(link = "log"),
REML = TRUE,
control = glmmTMBControl(
optimizer = optim,
optArgs = list(method = "BFGS")
)
)
AICc(copper.glmmTMB4a, copper.glmmTMB4b)
Conclusions:
copper.glmmTMB5a <- glmmTMB(WORMS ~ COPPER * DIST + (1 | PLATE),
dispformula = ~ COPPER * DIST,
data = copper,
family = tweedie(link = "log"),
REML = TRUE
)
copper.glmmTMB5b <- glmmTMB(WORMS ~ COPPER * DIST + (DIST | PLATE),
dispformula = ~ COPPER * DIST,
data = copper,
family = tweedie(link = "log"),
REML = TRUE
)
AICc(copper.glmmTMB5a, copper.glmmTMB5b)
Conclusions:
copper.glmmTMB6a <- glmmTMB(WORMS ~ COPPER * DIST + (1 | PLATE),
zi = ~1,
dispformula = ~ COPPER * DIST,
data = copper,
family = ziGamma(link = "log"),
REML = TRUE
)
copper.glmmTMB6b <- glmmTMB(WORMS ~ COPPER * DIST + (DIST | PLATE),
zi = ~1,
dispformula = ~ COPPER * DIST,
data = copper,
family = ziGamma(link = "log"),
REML = TRUE,
control = glmmTMBControl(
optimizer = optim,
optArgs = list(method = "BFGS")
)
)
AICc(copper.glmmTMB6a, copper.glmmTMB6b)
Conclusions:
plot_model(copper.glmmTMB1a, type = "diag")[-2] %>% plot_grid()
Conclusions:
copper.glmmTMB1a %>% performance::check_model()
Conclusions:
It is also possible to predict which modelling family would be the most suitable for the data. This can be attempted via an experimental routine in the performance package that uses random forests to classify a range of possible distributions.
copper.glmmTMB1a %>% performance::check_distribution()
Conclusions:
copper.resid <- copper.glmmTMB1a %>% simulateResiduals(plot = TRUE, integerResponse = TRUE)
Conclusions:
plot_model(copper.glmmTMB2a, type = "diag")[-2] %>% plot_grid()
Conclusions:
copper.glmmTMB2a %>% performance::check_model()
Conclusions:
It is also possible to predict which modelling family would be the most suitable for the data. This can be attempted via an experimental routine in the performance package that uses random forests to classify a range of possible distributions.
copper.glmmTMB2a %>% performance::check_distribution()
Conclusions:
copper.resid <- copper.glmmTMB2a %>% simulateResiduals(plot = TRUE, integerResponse = FALSE)
Conclusions:
plot_model(copper.glmmTMB3a, type = "diag")[-2] %>% plot_grid()
Conclusions:
copper.glmmTMB3a %>% performance::check_model()
Conclusions:
It is also possible to predict which modelling family would be the most suitable for the data. This can be attempted via an experimental routine in the performance package that uses random forests to classify a range of possible distributions.
copper.glmmTMB3a %>% performance::check_distribution()
Conclusions:
copper.resid <- copper.glmmTMB3a %>% simulateResiduals(plot = TRUE, integerResponse = FALSE)
Conclusions:
plot_model(copper.glmmTMB4a, type = "diag")[-2] %>% plot_grid()
Conclusions:
copper.glmmTMB4a %>% performance::check_model()
Conclusions:
It is also possible to predict which modelling family would be the most suitable for the data. This can be attempted via an experimental routine in the performance package that uses random forests to classify a range of possible distributions.
copper.glmmTMB4a %>% performance::check_distribution()
Conclusions:
copper.resid <- copper.glmmTMB4a %>% simulateResiduals(plot = TRUE, integerResponse = TRUE)
Conclusions:
plot_model(copper.glmmTMB5a, type = "diag")[-2] %>% plot_grid()
Conclusions:
copper.glmmTMB5a %>% performance::check_model()
Conclusions:
It is also possible to predict which modelling family would be the most suitable for the data. This can be attempted via an experimental routine in the performance package that uses random forests to classify a range of possible distributions.
copper.glmmTMB5a %>% performance::check_distribution()
Conclusions:
copper.resid <- copper.glmmTMB5a %>% simulateResiduals(plot = TRUE, integerResponse = TRUE)
Conclusions:
copper.resid <- copper.glmmTMB6a %>% simulateResiduals(plot = TRUE, integerResponse = TRUE)
Conclusions:
Conclusions:
copper.glmmTMB1a %>% plot_model(type = "eff", terms = c("DIST", "COPPER"))
Conclusions:
copper.glmmTMB1a %>%
allEffects() %>%
plot(multiline = TRUE, ci.style = "bars")
Conclusions:
copper.glmmTMB1a %>%
ggpredict(c("DIST", "COPPER")) %>%
plot()
Conclusions:
copper.glmmTMB1a %>%
ggemmeans(~ DIST + COPPER) %>%
plot()
Conclusions:
copper.glmmTMB2a %>% plot_model(type = "eff", terms = c("DIST", "COPPER"))
Conclusions:
Week 2 Dist 1 lower confidence interval extends below 0, if we were to back-transform the negatives, they will become positive and distort the intervalcopper.glmmTMB2a %>%
allEffects() %>%
plot(multiline = TRUE, ci.style = "bars")
Conclusions:
Week 2 Dist 1 lower confidence interval extends below 0, if we were to back-transform the negatives, they will become positive and distort the intervalcopper.glmmTMB2a %>%
ggpredict(c("DIST", "COPPER")) %>%
plot()
Conclusions:
Week 2 Dist 1 lower confidence interval have been distorted.copper.glmmTMB2a %>%
ggemmeans(~ DIST + COPPER) %>%
plot()
Conclusions:
Week 2 Dist 1 lower confidence interval have been distorted.plot_model(copper.glmmTMB3a, type = "eff", terms = c("DIST", "COPPER"))
copper.glmmTMB3a %>%
allEffects() %>%
plot(multiline = TRUE, ci.style = "bars")
copper.glmmTMB3a %>%
allEffects(transformation = NULL) %>%
plot(multiline = TRUE, ci.style = "bars")
copper.glmmTMB3a %>%
ggpredict(c("DIST", "COPPER")) %>%
plot()
copper.glmmTMB3a %>%
ggemmeans(~ DIST + COPPER) %>%
plot()
copper.glmmTMB4a %>% plot_model(type = "eff", terms = c("DIST", "COPPER"))
copper.glmmTMB4a %>%
allEffects() %>%
plot(multiline = TRUE, ci.style = "bars")
copper.glmmTMB4a %>%
allEffects(transformation = NULL) %>%
plot(multiline = TRUE, ci.style = "bars")
copper.glmmTMB4a %>%
ggpredict(c("DIST", "COPPER")) %>%
plot()
copper.glmmTMB4a %>%
ggemmeans(~ DIST + COPPER) %>%
plot()
copper.glmmTMB5a %>% plot_model(type = "eff", terms = c("DIST", "COPPER"))
copper.glmmTMB5a %>%
allEffects() %>%
plot(multiline = TRUE, ci.style = "bars")
copper.glmmTMB5a %>%
allEffects(transformation = NULL) %>%
plot(multiline = TRUE, ci.style = "bars")
copper.glmmTMB5a %>%
ggpredict(c("DIST", "COPPER")) %>%
plot()
copper.glmmTMB5a %>%
ggemmeans(~ DIST + COPPER) %>%
plot()
copper.glmmTMB1a %>% summary()
## Family: gaussian ( identity )
## Formula: WORMS ~ COPPER * DIST + (1 | PLATE)
## Dispersion: ~COPPER * DIST
## Data: copper
##
## AIC BIC logLik deviance df.resid
## 223.6 276.0 -86.8 173.6 47
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## PLATE (Intercept) 0.583 0.7636
## Residual NA NA
## Number of obs: 60, groups: PLATE, 15
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 10.8500 0.4041 26.851 < 2e-16 ***
## COPPERWeek 1 -3.6000 0.8008 -4.496 6.94e-06 ***
## COPPERWeek 2 -10.6000 0.5291 -20.036 < 2e-16 ***
## DIST2 1.1500 0.2920 3.938 8.22e-05 ***
## DIST3 1.5500 0.5789 2.677 0.007421 **
## DIST4 2.7000 1.2644 2.135 0.032729 *
## COPPERWeek 1:DIST2 -0.0500 1.0436 -0.048 0.961789
## COPPERWeek 2:DIST2 0.0500 0.4476 0.112 0.911061
## COPPERWeek 1:DIST3 -0.3000 1.0027 -0.299 0.764792
## COPPERWeek 2:DIST3 2.2000 0.7851 2.802 0.005078 **
## COPPERWeek 1:DIST4 0.0500 1.5472 0.032 0.974219
## COPPERWeek 2:DIST4 4.9000 1.3656 3.588 0.000333 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.4550 1.4707 -0.989 0.3225
## COPPERWeek 1 2.0465 1.9052 1.074 0.2827
## COPPERWeek 2 -9.5217 87.7945 -0.108 0.9136
## DIST2 -0.1899 1.4619 -0.130 0.8967
## DIST3 1.8213 1.8955 0.961 0.3366
## DIST4 3.5040 1.7005 2.061 0.0393 *
## COPPERWeek 1:DIST2 0.7655 1.8330 0.418 0.6762
## COPPERWeek 2:DIST2 10.6139 87.7926 0.121 0.9038
## COPPERWeek 1:DIST3 -1.9782 2.8271 -0.700 0.4841
## COPPERWeek 2:DIST3 9.4963 87.8062 0.108 0.9139
## COPPERWeek 1:DIST4 -3.3216 2.6021 -1.276 0.2018
## COPPERWeek 2:DIST4 7.7583 87.8018 0.088 0.9296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conclusions:
copper.glmmTMB1a %>% tidy(conf.int = TRUE)
copper.glmmTMB1a %>%
tidy(conf.int = TRUE) %>%
kable()
| effect | component | group | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 10.8500000 | 0.4040825 | 26.8509520 | 0.0000000 | 10.0580128 | 11.6419872 |
| fixed | cond | NA | COPPERWeek 1 | -3.6000000 | 0.8007748 | -4.4956459 | 0.0000069 | -5.1694898 | -2.0305102 |
| fixed | cond | NA | COPPERWeek 2 | -10.6000000 | 0.5290450 | -20.0361014 | 0.0000000 | -11.6369092 | -9.5630908 |
| fixed | cond | NA | DIST2 | 1.1500000 | 0.2920393 | 3.9378257 | 0.0000822 | 0.5776134 | 1.7223866 |
| fixed | cond | NA | DIST3 | 1.5500000 | 0.5789338 | 2.6773353 | 0.0074210 | 0.4153105 | 2.6846895 |
| fixed | cond | NA | DIST4 | 2.7000000 | 1.2644060 | 2.1353900 | 0.0327292 | 0.2218097 | 5.1781903 |
| fixed | cond | NA | COPPERWeek 1:DIST2 | -0.0500000 | 1.0436511 | -0.0479087 | 0.9617890 | -2.0955187 | 1.9955187 |
| fixed | cond | NA | COPPERWeek 2:DIST2 | 0.0500000 | 0.4476264 | 0.1117003 | 0.9110610 | -0.8273316 | 0.9273316 |
| fixed | cond | NA | COPPERWeek 1:DIST3 | -0.3000000 | 1.0026934 | -0.2991941 | 0.7647919 | -2.2652430 | 1.6652430 |
| fixed | cond | NA | COPPERWeek 2:DIST3 | 2.2000000 | 0.7851341 | 2.8020691 | 0.0050776 | 0.6611654 | 3.7388346 |
| fixed | cond | NA | COPPERWeek 1:DIST4 | 0.0500000 | 1.5471737 | 0.0323170 | 0.9742193 | -2.9824048 | 3.0824048 |
| fixed | cond | NA | COPPERWeek 2:DIST4 | 4.9000000 | 1.3655898 | 3.5881933 | 0.0003330 | 2.2234932 | 7.5765068 |
| ran_pars | cond | PLATE | sd__(Intercept) | 0.7635527 | NA | NA | NA | 0.4030681 | 1.4464373 |
| ran_pars | cond | Residual | sd__Observation | NA | NA | NA | NA | 0.4030681 | 1.4464373 |
Conclusions:
# warning this is only appropriate for html output
copper.glmmTMB1a %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| Â | WORMS | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 10.85 | 0.40 | 10.06 – 11.64 | <0.001 |
| COPPER [Week 1] | -3.60 | 0.80 | -5.17 – -2.03 | <0.001 |
| COPPER [Week 2] | -10.60 | 0.53 | -11.64 – -9.56 | <0.001 |
| DIST [2] | 1.15 | 0.29 | 0.58 – 1.72 | <0.001 |
| DIST [3] | 1.55 | 0.58 | 0.42 – 2.68 | 0.007 |
| DIST [4] | 2.70 | 1.26 | 0.22 – 5.18 | 0.033 |
|
COPPER [Week 1] * DIST [2] |
-0.05 | 1.04 | -2.10 – 2.00 | 0.962 |
|
COPPER [Week 2] * DIST [2] |
0.05 | 0.45 | -0.83 – 0.93 | 0.911 |
|
COPPER [Week 1] * DIST [3] |
-0.30 | 1.00 | -2.27 – 1.67 | 0.765 |
|
COPPER [Week 2] * DIST [3] |
2.20 | 0.79 | 0.66 – 3.74 | 0.005 |
|
COPPER [Week 1] * DIST [4] |
0.05 | 1.55 | -2.98 – 3.08 | 0.974 |
|
COPPER [Week 2] * DIST [4] |
4.90 | 1.37 | 2.22 – 7.58 | <0.001 |
| Random Effects | ||||
| σ2 | NA | |||
| τ00 PLATE | 0.58 | |||
| N PLATE | 15 | |||
| Observations | 60 | |||
| AIC | 223.648 | |||
Conclusions:
copper.glmmTMB2a %>% summary()
## Family: gaussian ( identity )
## Formula: sqrt(WORMS) ~ COPPER * DIST + (1 | PLATE)
## Dispersion: ~COPPER * DIST
## Data: copper
##
## AIC BIC logLik deviance df.resid
## 70.4 122.8 -10.2 20.4 47
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## PLATE (Intercept) 0.03425 0.1851
## Residual NA NA
## Number of obs: 60, groups: PLATE, 15
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.288733 0.086711 37.93 < 2e-16 ***
## COPPERWeek 1 -0.603217 0.205378 -2.94 0.00331 **
## COPPERWeek 2 -2.988733 0.220762 -13.54 < 2e-16 ***
## DIST2 0.171972 0.036295 4.74 2.16e-06 ***
## DIST3 0.229034 0.088727 2.58 0.00984 **
## DIST4 0.377765 0.173294 2.18 0.02926 *
## COPPERWeek 1:DIST2 0.018871 0.251141 0.08 0.94010
## COPPERWeek 2:DIST2 0.671908 0.282941 2.37 0.01756 *
## COPPERWeek 1:DIST3 -0.007477 0.190562 -0.04 0.96870
## COPPERWeek 2:DIST3 1.456036 0.228078 6.38 1.73e-10 ***
## COPPERWeek 1:DIST4 0.091074 0.243731 0.37 0.70865
## COPPERWeek 2:DIST4 2.119555 0.264172 8.02 1.03e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.70185 1.44255 -3.953 7.73e-05 ***
## COPPERWeek 1 3.72894 1.62471 2.295 0.02173 *
## COPPERWeek 2 3.94060 1.63297 2.413 0.01582 *
## DIST2 -0.02827 2.49264 -0.011 0.99095
## DIST3 2.37823 1.60298 1.484 0.13791
## DIST4 3.78325 1.52792 2.476 0.01328 *
## COPPERWeek 1:DIST2 0.22759 2.69205 0.085 0.93263
## COPPERWeek 2:DIST2 0.28382 2.68454 0.106 0.91580
## COPPERWeek 1:DIST3 -6.16395 2.92030 -2.111 0.03480 *
## COPPERWeek 2:DIST3 -3.63483 1.98292 -1.833 0.06679 .
## COPPERWeek 1:DIST4 -6.66142 1.99284 -3.343 0.00083 ***
## COPPERWeek 2:DIST4 -5.63576 2.19170 -2.571 0.01013 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conclusions:
copper.glmmTMB2a %>% tidy(conf.int = TRUE)
copper.glmmTMB2a %>%
tidy(conf.int = TRUE) %>%
kable()
| effect | component | group | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 3.2887330 | 0.0867111 | 37.9274768 | 0.0000000 | 3.1187823 | 3.4586836 |
| fixed | cond | NA | COPPERWeek 1 | -0.6032170 | 0.2053778 | -2.9371094 | 0.0033129 | -1.0057500 | -0.2006839 |
| fixed | cond | NA | COPPERWeek 2 | -2.9887330 | 0.2207616 | -13.5382851 | 0.0000000 | -3.4214177 | -2.5560482 |
| fixed | cond | NA | DIST2 | 0.1719716 | 0.0362946 | 4.7382140 | 0.0000022 | 0.1008355 | 0.2431077 |
| fixed | cond | NA | DIST3 | 0.2290344 | 0.0887268 | 2.5813438 | 0.0098417 | 0.0551330 | 0.4029357 |
| fixed | cond | NA | DIST4 | 0.3777650 | 0.1732937 | 2.1799124 | 0.0292640 | 0.0381157 | 0.7174144 |
| fixed | cond | NA | COPPERWeek 1:DIST2 | 0.0188707 | 0.2511415 | 0.0751395 | 0.9401037 | -0.4733576 | 0.5110989 |
| fixed | cond | NA | COPPERWeek 2:DIST2 | 0.6719077 | 0.2829410 | 2.3747274 | 0.0175619 | 0.1173536 | 1.2264619 |
| fixed | cond | NA | COPPERWeek 1:DIST3 | -0.0074774 | 0.1905620 | -0.0392389 | 0.9686999 | -0.3809722 | 0.3660173 |
| fixed | cond | NA | COPPERWeek 2:DIST3 | 1.4560361 | 0.2280783 | 6.3839301 | 0.0000000 | 1.0090108 | 1.9030615 |
| fixed | cond | NA | COPPERWeek 1:DIST4 | 0.0910738 | 0.2437315 | 0.3736643 | 0.7086541 | -0.3866312 | 0.5687787 |
| fixed | cond | NA | COPPERWeek 2:DIST4 | 2.1195548 | 0.2641719 | 8.0233926 | 0.0000000 | 1.6017874 | 2.6373222 |
| ran_pars | cond | PLATE | sd__(Intercept) | 0.1850791 | NA | NA | NA | 0.1144146 | 0.2993873 |
| ran_pars | cond | Residual | sd__Observation | NA | NA | NA | NA | 0.1144146 | 0.2993873 |
Conclusions:
# warning this is only appropriate for html output
copper.glmmTMB2a %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| Â | sqrt(WORMS) | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 3.29 | 0.09 | 3.12 – 3.46 | <0.001 |
| COPPER [Week 1] | -0.60 | 0.21 | -1.01 – -0.20 | 0.003 |
| COPPER [Week 2] | -2.99 | 0.22 | -3.42 – -2.56 | <0.001 |
| DIST [2] | 0.17 | 0.04 | 0.10 – 0.24 | <0.001 |
| DIST [3] | 0.23 | 0.09 | 0.06 – 0.40 | 0.010 |
| DIST [4] | 0.38 | 0.17 | 0.04 – 0.72 | 0.029 |
|
COPPER [Week 1] * DIST [2] |
0.02 | 0.25 | -0.47 – 0.51 | 0.940 |
|
COPPER [Week 2] * DIST [2] |
0.67 | 0.28 | 0.12 – 1.23 | 0.018 |
|
COPPER [Week 1] * DIST [3] |
-0.01 | 0.19 | -0.38 – 0.37 | 0.969 |
|
COPPER [Week 2] * DIST [3] |
1.46 | 0.23 | 1.01 – 1.90 | <0.001 |
|
COPPER [Week 1] * DIST [4] |
0.09 | 0.24 | -0.39 – 0.57 | 0.709 |
|
COPPER [Week 2] * DIST [4] |
2.12 | 0.26 | 1.60 – 2.64 | <0.001 |
| Random Effects | ||||
| σ2 | NA | |||
| τ00 PLATE | 0.03 | |||
| N PLATE | 15 | |||
| Observations | 60 | |||
| AIC | 70.420 | |||
Conclusions:
copper.glmmTMB3a %>% summary()
## Family: gaussian ( log )
## Formula: WORMSp ~ COPPER * DIST + (1 | PLATE)
## Dispersion: ~COPPER * DIST
## Data: copper
##
## AIC BIC logLik deviance df.resid
## 261.8 314.2 -105.9 211.8 47
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## PLATE (Intercept) 0.01481 0.1217
## Residual NA NA
## Number of obs: 60, groups: PLATE, 15
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.38213 0.05716 41.67 < 2e-16 ***
## COPPERWeek 1 -0.43674 0.14703 -2.97 0.00297 **
## COPPERWeek 2 -3.53420 0.55522 -6.37 1.95e-10 ***
## DIST2 0.09842 0.02367 4.16 3.21e-05 ***
## DIST3 0.12737 0.05013 2.54 0.01106 *
## DIST4 0.21529 0.09627 2.24 0.02533 *
## COPPERWeek 1:DIST2 0.04587 0.18276 0.25 0.80181
## COPPERWeek 2:DIST2 1.38039 0.65684 2.10 0.03559 *
## COPPERWeek 1:DIST3 0.06134 0.13507 0.45 0.64973
## COPPERWeek 2:DIST3 2.41188 0.56303 4.28 1.84e-05 ***
## COPPERWeek 1:DIST4 0.13438 0.15817 0.85 0.39557
## COPPERWeek 2:DIST4 2.99625 0.55827 5.37 8.01e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.72906 1.33324 -1.297 0.19467
## COPPERWeek 1 3.08496 1.52096 2.028 0.04253 *
## COPPERWeek 2 -0.15852 1.51377 -0.105 0.91660
## DIST2 0.06408 2.38375 0.027 0.97855
## DIST3 2.26009 1.54021 1.467 0.14227
## DIST4 3.83419 1.44032 2.662 0.00777 **
## COPPERWeek 1:DIST2 0.35915 2.58817 0.139 0.88963
## COPPERWeek 2:DIST2 2.06889 2.58080 0.802 0.42276
## COPPERWeek 1:DIST3 -5.54247 2.23181 -2.483 0.01301 *
## COPPERWeek 2:DIST3 -0.28135 1.85305 -0.152 0.87932
## COPPERWeek 1:DIST4 -6.75562 2.13996 -3.157 0.00159 **
## COPPERWeek 2:DIST4 -3.16616 3.15813 -1.002 0.31608
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conclusions:
copper.glmmTMB3a %>% tidy(conf.int = TRUE)
## or on the response scale
copper.glmmTMB3a %>% tidy(conf.int = TRUE, exponentiate = TRUE)
copper.glmmTMB3a %>%
tidy(conf.int = TRUE, exponentiate = TRUE) %>%
kable()
| effect | component | group | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 10.8279548 | 0.6189666 | 41.6720531 | 0.0000000 | 9.6802941 | 12.1116779 |
| fixed | cond | NA | COPPERWeek 1 | 0.6461422 | 0.0950040 | -2.9703326 | 0.0029748 | 0.4843661 | 0.8619508 |
| fixed | cond | NA | COPPERWeek 2 | 0.0291820 | 0.0162023 | -6.3654628 | 0.0000000 | 0.0098291 | 0.0866394 |
| fixed | cond | NA | DIST2 | 1.1034259 | 0.0261172 | 4.1581387 | 0.0000321 | 1.0534063 | 1.1558206 |
| fixed | cond | NA | DIST3 | 1.1358411 | 0.0569415 | 2.5407803 | 0.0110605 | 1.0295453 | 1.2531114 |
| fixed | cond | NA | DIST4 | 1.2402272 | 0.1194020 | 2.2362620 | 0.0253346 | 1.0269573 | 1.4977872 |
| fixed | cond | NA | COPPERWeek 1:DIST2 | 1.0469424 | 0.1913425 | 0.2510018 | 0.8018127 | 0.7317361 | 1.4979286 |
| fixed | cond | NA | COPPERWeek 2:DIST2 | 3.9764414 | 2.6118909 | 2.1015538 | 0.0355924 | 1.0974640 | 14.4078410 |
| fixed | cond | NA | COPPERWeek 1:DIST3 | 1.0632593 | 0.1436099 | 0.4541418 | 0.6497268 | 0.8159640 | 1.3855027 |
| fixed | cond | NA | COPPERWeek 2:DIST3 | 11.1549560 | 6.2805624 | 4.2837658 | 0.0000184 | 3.7001285 | 33.6293844 |
| fixed | cond | NA | COPPERWeek 1:DIST4 | 1.1438246 | 0.1809206 | 0.8495681 | 0.3955652 | 0.8389257 | 1.5595357 |
| fixed | cond | NA | COPPERWeek 2:DIST4 | 20.0103535 | 11.1712358 | 5.3669996 | 0.0000001 | 6.6996499 | 59.7664441 |
| ran_pars | cond | PLATE | sd__(Intercept) | 0.1217129 | NA | NA | NA | 0.0750934 | 0.1972747 |
| ran_pars | cond | Residual | sd__Observation | NA | NA | NA | NA | 0.0750934 | 0.1972747 |
Conclusions:
# warning this is only appropriate for html output
copper.glmmTMB3a %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| Â | WORMSp | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 2.38 | 0.06 | 2.27 – 2.49 | <0.001 |
| COPPER [Week 1] | -0.44 | 0.15 | -0.72 – -0.15 | 0.003 |
| COPPER [Week 2] | -3.53 | 0.56 | -4.62 – -2.45 | <0.001 |
| DIST [2] | 0.10 | 0.02 | 0.05 – 0.14 | <0.001 |
| DIST [3] | 0.13 | 0.05 | 0.03 – 0.23 | 0.011 |
| DIST [4] | 0.22 | 0.10 | 0.03 – 0.40 | 0.025 |
|
COPPER [Week 1] * DIST [2] |
0.05 | 0.18 | -0.31 – 0.40 | 0.802 |
|
COPPER [Week 2] * DIST [2] |
1.38 | 0.66 | 0.09 – 2.67 | 0.036 |
|
COPPER [Week 1] * DIST [3] |
0.06 | 0.14 | -0.20 – 0.33 | 0.650 |
|
COPPER [Week 2] * DIST [3] |
2.41 | 0.56 | 1.31 – 3.52 | <0.001 |
|
COPPER [Week 1] * DIST [4] |
0.13 | 0.16 | -0.18 – 0.44 | 0.396 |
|
COPPER [Week 2] * DIST [4] |
3.00 | 0.56 | 1.90 – 4.09 | <0.001 |
| Random Effects | ||||
| σ2 | NA | |||
| τ00 PLATE | 0.01 | |||
| N PLATE | 15 | |||
| Observations | 60 | |||
| AIC | 261.825 | |||
copper.glmmTMB4a %>% summary()
## Family: Gamma ( log )
## Formula: WORMSp ~ COPPER * DIST + (1 | PLATE)
## Dispersion: ~COPPER * DIST
## Data: copper
##
## AIC BIC logLik deviance df.resid
## 254.8 307.1 -102.4 204.8 47
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## PLATE (Intercept) 0.01362 0.1167
## Number of obs: 60, groups: PLATE, 15
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.378243 0.055023 43.22 < 2e-16 ***
## COPPERWeek 1 -0.382097 0.136433 -2.80 0.00510 **
## COPPERWeek 2 -3.492196 0.428225 -8.16 3.49e-16 ***
## DIST2 0.102890 0.022118 4.65 3.29e-06 ***
## DIST3 0.140327 0.053928 2.60 0.00926 **
## DIST4 0.230374 0.096130 2.40 0.01655 *
## COPPERWeek 1:DIST2 0.034258 0.167542 0.20 0.83798
## COPPERWeek 2:DIST2 1.413100 0.534279 2.64 0.00817 **
## COPPERWeek 1:DIST3 -0.007238 0.127867 -0.06 0.95486
## COPPERWeek 2:DIST3 2.351935 0.437734 5.37 7.74e-08 ***
## COPPERWeek 1:DIST4 0.066969 0.151014 0.44 0.65743
## COPPERWeek 2:DIST4 2.938791 0.434155 6.77 1.30e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.48905 1.21527 5.340 9.32e-08 ***
## COPPERWeek 1 -3.74194 1.43230 -2.613 0.008987 **
## COPPERWeek 2 -6.36019 1.36688 -4.653 3.27e-06 ***
## DIST2 0.50194 2.62420 0.191 0.848310
## DIST3 -2.14711 1.40115 -1.532 0.125425
## DIST4 -3.37965 1.33539 -2.531 0.011380 *
## COPPERWeek 1:DIST2 -0.64196 2.80996 -0.228 0.819290
## COPPERWeek 2:DIST2 -0.02106 2.76701 -0.008 0.993928
## COPPERWeek 1:DIST3 5.25749 2.24455 2.342 0.019163 *
## COPPERWeek 2:DIST3 4.90007 1.71157 2.863 0.004198 **
## COPPERWeek 1:DIST4 6.30087 2.07341 3.039 0.002375 **
## COPPERWeek 2:DIST4 8.19547 2.44890 3.347 0.000818 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conclusions:
copper.glmmTMB4a %>% tidy(conf.int = TRUE)
## or on the response scale
copper.glmmTMB4a %>% tidy(conf.int = TRUE, exponentiate = TRUE)
copper.glmmTMB4a %>%
tidy(conf.int = TRUE, exponentiate = TRUE) %>%
kable()
| effect | component | group | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 10.7859394 | 0.5934698 | 43.2230751 | 0.0000000 | 9.6832848 | 12.0141554 |
| fixed | cond | NA | COPPERWeek 1 | 0.6824286 | 0.0931055 | -2.8006302 | 0.0051003 | 0.5223067 | 0.8916387 |
| fixed | cond | NA | COPPERWeek 2 | 0.0304340 | 0.0130326 | -8.1550569 | 0.0000000 | 0.0131478 | 0.0704473 |
| fixed | cond | NA | DIST2 | 1.1083698 | 0.0245154 | 4.6517970 | 0.0000033 | 1.0613472 | 1.1574757 |
| fixed | cond | NA | DIST3 | 1.1506504 | 0.0620523 | 2.6021240 | 0.0092648 | 1.0352370 | 1.2789307 |
| fixed | cond | NA | DIST4 | 1.2590708 | 0.1210345 | 2.3964843 | 0.0165532 | 1.0428556 | 1.5201140 |
| fixed | cond | NA | COPPERWeek 1:DIST2 | 1.0348512 | 0.1733810 | 0.2044722 | 0.8379845 | 0.7451882 | 1.4371095 |
| fixed | cond | NA | COPPERWeek 2:DIST2 | 4.1086720 | 2.1951767 | 2.6448731 | 0.0081722 | 1.4418579 | 11.7079401 |
| fixed | cond | NA | COPPERWeek 1:DIST3 | 0.9927886 | 0.1269450 | -0.0566023 | 0.9548620 | 0.7727092 | 1.2755499 |
| fixed | cond | NA | COPPERWeek 2:DIST3 | 10.5058839 | 4.5987817 | 5.3729797 | 0.0000001 | 4.4548472 | 24.7760678 |
| fixed | cond | NA | COPPERWeek 1:DIST4 | 1.0692628 | 0.1614733 | 0.4434661 | 0.6574286 | 0.7953190 | 1.4375652 |
| fixed | cond | NA | COPPERWeek 2:DIST4 | 18.8929929 | 8.2024873 | 6.7689906 | 0.0000000 | 8.0676557 | 44.2439780 |
| ran_pars | cond | PLATE | sd__(Intercept) | 0.1166844 | NA | NA | NA | 0.0699492 | 0.1946447 |
Conclusions:
# warning this is only appropriate for html output
copper.glmmTMB4a %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| Â | WORMSp | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 10.79 | 0.59 | 9.68 – 12.01 | <0.001 |
| COPPER [Week 1] | 0.68 | 0.09 | 0.52 – 0.89 | 0.005 |
| COPPER [Week 2] | 0.03 | 0.01 | 0.01 – 0.07 | <0.001 |
| DIST [2] | 1.11 | 0.02 | 1.06 – 1.16 | <0.001 |
| DIST [3] | 1.15 | 0.06 | 1.04 – 1.28 | 0.009 |
| DIST [4] | 1.26 | 0.12 | 1.04 – 1.52 | 0.017 |
|
COPPER [Week 1] * DIST [2] |
1.03 | 0.17 | 0.75 – 1.44 | 0.838 |
|
COPPER [Week 2] * DIST [2] |
4.11 | 2.20 | 1.44 – 11.71 | 0.008 |
|
COPPER [Week 1] * DIST [3] |
0.99 | 0.13 | 0.77 – 1.28 | 0.955 |
|
COPPER [Week 2] * DIST [3] |
10.51 | 4.60 | 4.45 – 24.78 | <0.001 |
|
COPPER [Week 1] * DIST [4] |
1.07 | 0.16 | 0.80 – 1.44 | 0.657 |
|
COPPER [Week 2] * DIST [4] |
18.89 | 8.20 | 8.07 – 44.24 | <0.001 |
| Random Effects | ||||
| σ2 | NA | |||
| τ00 PLATE | 0.01 | |||
| N PLATE | 15 | |||
| Observations | 60 | |||
| AIC | 254.770 | |||
copper.glmmTMB5a %>% summary()
## Family: tweedie ( log )
## Formula: WORMS ~ COPPER * DIST + (1 | PLATE)
## Dispersion: ~COPPER * DIST
## Data: copper
##
## AIC BIC logLik deviance df.resid
## 266.2 320.7 -107.1 214.2 46
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## PLATE (Intercept) 0.014 0.1183
## Number of obs: 60, groups: PLATE, 15
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.37938 0.05573 42.70 < 2e-16 ***
## COPPERWeek 1 -0.39612 0.13909 -2.85 0.004399 **
## COPPERWeek 2 -3.76158 0.88806 -4.24 2.28e-05 ***
## DIST2 0.10166 0.02263 4.49 7.03e-06 ***
## DIST3 0.13629 0.05283 2.58 0.009881 **
## DIST4 0.22547 0.09579 2.35 0.018581 *
## COPPERWeek 1:DIST2 0.03799 0.17062 0.22 0.823814
## COPPERWeek 2:DIST2 1.66468 0.93924 1.77 0.076333 .
## COPPERWeek 1:DIST3 0.01145 0.12965 0.09 0.929646
## COPPERWeek 2:DIST3 2.62722 0.89294 2.94 0.003259 **
## COPPERWeek 1:DIST4 0.08549 0.15257 0.56 0.575260
## COPPERWeek 2:DIST4 3.21325 0.89069 3.61 0.000309 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Dispersion model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.1002 1.5636 -3.262 0.00111 **
## COPPERWeek 1 3.5513 1.4528 2.444 0.01451 *
## COPPERWeek 2 5.6591 2.6184 2.161 0.03067 *
## DIST2 -0.3254 2.5106 -0.130 0.89687
## DIST3 2.1778 1.4358 1.517 0.12932
## DIST4 3.5035 1.3621 2.572 0.01011 *
## COPPERWeek 1:DIST2 0.5416 2.6998 0.201 0.84102
## COPPERWeek 2:DIST2 -0.7091 2.9326 -0.242 0.80893
## COPPERWeek 1:DIST3 -5.3508 2.2101 -2.421 0.01547 *
## COPPERWeek 2:DIST3 -4.7669 2.5198 -1.892 0.05852 .
## COPPERWeek 1:DIST4 -6.4231 2.0580 -3.121 0.00180 **
## COPPERWeek 2:DIST4 -7.9456 3.3067 -2.403 0.01627 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conclusions:
copper.glmmTMB5a %>% tidy(conf.int = FALSE)
## or on the response scale
copper.glmmTMB5a %>% tidy(conf.int = FALSE, exponentiate = TRUE)
copper.glmmTMB5a %>%
tidy(conf.int = FALSE, exponentiate = TRUE) %>%
kable()
| effect | component | group | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 10.7981605 | 0.6017546 | 42.6966095 | 0.0000000 |
| fixed | cond | NA | COPPERWeek 1 | 0.6729242 | 0.0935940 | -2.8480520 | 0.0043988 |
| fixed | cond | NA | COPPERWeek 2 | 0.0232470 | 0.0206448 | -4.2357183 | 0.0000228 |
| fixed | cond | NA | DIST2 | 1.1070069 | 0.0250480 | 4.4928967 | 0.0000070 |
| fixed | cond | NA | DIST3 | 1.1460102 | 0.0605381 | 2.5799554 | 0.0098813 |
| fixed | cond | NA | DIST4 | 1.2529125 | 0.1200148 | 2.3538379 | 0.0185807 |
| fixed | cond | NA | COPPERWeek 1:DIST2 | 1.0387185 | 0.1772289 | 0.2226419 | 0.8238142 |
| fixed | cond | NA | COPPERWeek 2:DIST2 | 5.2839933 | 4.9629415 | 1.7723701 | 0.0763331 |
| fixed | cond | NA | COPPERWeek 1:DIST3 | 1.0115129 | 0.1311473 | 0.0882897 | 0.9296464 |
| fixed | cond | NA | COPPERWeek 2:DIST3 | 13.8352323 | 12.3540811 | 2.9422000 | 0.0032589 |
| fixed | cond | NA | COPPERWeek 1:DIST4 | 1.0892460 | 0.1661812 | 0.5603216 | 0.5752601 |
| fixed | cond | NA | COPPERWeek 2:DIST4 | 24.8597060 | 22.1422778 | 3.6075966 | 0.0003090 |
| ran_pars | cond | PLATE | sd__(Intercept) | 0.1183171 | NA | NA | NA |
Conclusions:
# warning this is only appropriate for html output
copper.glmmTMB5a %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| Â | WORMS | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 10.80 | 0.60 | 9.68 – 12.04 | <0.001 |
| COPPER [Week 1] | 0.67 | 0.09 | 0.51 – 0.88 | 0.004 |
| COPPER [Week 2] | 0.02 | 0.02 | 0.00 – 0.13 | <0.001 |
| DIST [2] | 1.11 | 0.03 | 1.06 – 1.16 | <0.001 |
| DIST [3] | 1.15 | 0.06 | 1.03 – 1.27 | 0.010 |
| DIST [4] | 1.25 | 0.12 | 1.04 – 1.51 | 0.019 |
|
COPPER [Week 1] * DIST [2] |
1.04 | 0.18 | 0.74 – 1.45 | 0.824 |
|
COPPER [Week 2] * DIST [2] |
5.28 | 4.96 | 0.84 – 33.30 | 0.076 |
|
COPPER [Week 1] * DIST [3] |
1.01 | 0.13 | 0.78 – 1.30 | 0.930 |
|
COPPER [Week 2] * DIST [3] |
13.84 | 12.35 | 2.40 – 79.63 | 0.003 |
|
COPPER [Week 1] * DIST [4] |
1.09 | 0.17 | 0.81 – 1.47 | 0.575 |
|
COPPER [Week 2] * DIST [4] |
24.86 | 22.14 | 4.34 – 142.45 | <0.001 |
| Random Effects | ||||
| σ2 | NA | |||
| τ00 PLATE | 0.01 | |||
| N PLATE | 15 | |||
| Observations | 60 | |||
| AIC | 266.231 | |||
The presence of an interaction suggests that the effect of the copper treatment is not consistent across the distances and vice verse. We may wish to further explore the differences in copper treatments separately for each distance.
copper.glmmTMB4a %>%
emmeans(~ COPPER | DIST, type = "response") %>%
pairs()
## DIST = 1:
## contrast ratio SE df null t.ratio p.value
## control / Week 1 1.47 0.200 47 1 2.801 0.0198
## control / Week 2 32.86 14.071 47 1 8.155 <.0001
## Week 1 / Week 2 22.42 9.925 47 1 7.027 <.0001
##
## DIST = 2:
## contrast ratio SE df null t.ratio p.value
## control / Week 1 1.42 0.202 47 1 2.434 0.0484
## control / Week 2 8.00 2.720 47 1 6.113 <.0001
## Week 1 / Week 2 5.65 2.038 47 1 4.797 <.0001
##
## DIST = 3:
## contrast ratio SE df null t.ratio p.value
## control / Week 1 1.48 0.137 47 1 4.193 0.0004
## control / Week 2 3.13 0.435 47 1 8.206 <.0001
## Week 1 / Week 2 2.12 0.279 47 1 5.713 <.0001
##
## DIST = 4:
## contrast ratio SE df null t.ratio p.value
## control / Week 1 1.37 0.168 47 1 2.566 0.0354
## control / Week 2 1.74 0.219 47 1 4.400 0.0002
## Week 1 / Week 2 1.27 0.110 47 1 2.738 0.0232
##
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log scale
## copper.glmmTMB4a %>% emmeans(~COPPER|DIST, type='response') %>% contrast(method='pairwise')
Conclusions:
copper.glmmTMB4a %>%
emmeans(~ DIST | COPPER, type = "response") %>%
contrast(method = "poly")
## COPPER = control:
## contrast ratio SE df null t.ratio p.value
## linear 2.07 0.61 47 1 2.484 0.0166
## quadratic 0.99 0.11 47 1 -0.117 0.9073
## cubic 1.13 0.21 47 1 0.637 0.5270
##
## COPPER = Week 1:
## contrast ratio SE df null t.ratio p.value
## linear 2.43 0.90 47 1 2.393 0.0208
## quadratic 1.03 0.17 47 1 0.160 0.8739
## cubic 1.36 0.53 47 1 0.795 0.4306
##
## COPPER = Week 2:
## contrast ratio SE df null t.ratio p.value
## linear 35730.79 47166.69 47 1 7.942 <.0001
## quadratic 0.43 0.24 47 1 -1.538 0.1308
## cubic 1.27 1.43 47 1 0.214 0.8315
##
## Tests are performed on the log scale
## copper.glmmTMB4a %>% emmeans(~DIST|COPPER) %>% regrid() %>% contrast(method='poly')
Conclusions:
newdata <- copper.glmmTMB4a %>%
emmeans(~ COPPER | DIST, type = "response") %>%
summary(infer = TRUE) %>%
as.data.frame()
head(newdata)
g1 <- ggplot(newdata, aes(y = response, x = DIST, fill = COPPER)) +
geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL),
shape = 21,
position = position_dodge(width = 0.2)
) +
theme_bw()
copper.comp <- copper.glmmTMB4a %>%
emmeans(~ COPPER | DIST, type = "response") %>%
pairs() %>%
## contrast(method='pairwise') %>%
summary(infer = TRUE) %>%
as.data.frame()
head(copper.comp)
g2 <- ggplot(copper.comp, aes(y = ratio, x = contrast, color = DIST)) +
geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL),
position = position_dodge(width = 0.6)
) +
geom_hline(yintercept = 1) +
scale_y_continuous(trans = scales::log2_trans(), limits = c(0.25, 500)) +
coord_flip() +
theme_bw()
g1 + g2